clc;
clear;
format short

%% 网格离散数据
JM = [1 2; 2 3; 3 4; 4 5; 5 6]; % 单元数据
JX = [0; 0.2; 0.4; 0.6; 0.8; 1]; % 节点数据
JB1 = [1, 0]; % Dirichlet边界(节点1的边界数值为0)
N = 6; % 节点总数
E = 5; % 单元总数

%% 单元方程计算
ke = [-1/2 1/2; -1/2 1/2]; % 单元刚度矩阵
delta_x = 1/5; % 单元长度
be = delta_x * [1/2; 1/2]; % 单元荷载向量

%% 总体方程组合
K = zeros(6, 6); % 初始化总刚度矩阵 K, 节点总数为6
b = zeros(6, 1); % 初始化总载荷向量 b
for i = 1:E
    for m = 1:2 % 单元内双层循环
        for n = 1:2
            K(JM(i, m), JM(i, n)) = K(JM(i, m), JM(i, n)) + ke(m, n);
        end
    end
    for m = 1:2 % 单元内单层循环
        b(JM(i, m), 1) = b(JM(i, m), 1) + be(m, 1);
    end
end

%% 代入 JB1 求解
for i = 1:length(JB1(:, 1))
    II = JB1(i, 1); % 提取边界节点编号
    b = b - K(:, II)*JB1(i, 2); % 对角线归一法第一步
    K(II, :) = K(II, :)*0; % 对角线归一法第二步
    K(:, II) = K(:, II)*0; % 对角线归一法第三步
    K(II, II) = 1; % 对角线归一法第四步
    b(II, 1) = 0; % 对角线归一法第五步
end

%% 求解方程
u = K\b; % 求解方程
result = [JX, u] % 输出结果

%% 绘图
plot(JX, u, '*');

